/*----------------------------------------------------------------------------

Code for: 

	- Table 3
	- Figures 1, 2, 3, 4.

----------------------------------------------------------------------------*/


global replica_dir="<add path to your working directory>"

global data_dir=	"$replica_dir/data"
global output_dir=	"$replica_dir/output"

use "$data_dir/slvl_fin.dta", clear  



/*----------------------------------------------------------------------------*/

**# Fig 4. Density of evaluation index and homogeneous group cutoffs

/*----------------------------------------------------------------------------*/	


twoway 	(hist indicer,  w(0.2) fcolor(black) lcolor(black) ) ///
		(hist threshold_1_hg,  w(0.2) fcolor(gray) lcolor(gray) ), ///
		graphregion(color(white)) ///
		legend(label(1 "SNED Index") label(2 "100% bonus cutoff points") position(6)  region(lstyle(none))) ///
		ylabel(0(0.05)0.2)
		
graph export "$output_dir/hist_ind_cutoffs_25.png", as(png) replace 

twoway 	(hist indicer,  w(0.2) fcolor(black) lcolor(black) ) ///
		(hist threshold_2_hg,  w(0.2) fcolor(gray) lcolor(gray) ), ///
		graphregion(color(white)) ///
		legend(label(1 "SNED Index") label(2 "60% bonus cutoff points") position(6)  region(lstyle(none)))  ///
		ylabel(0(0.05)0.2) 
		
graph export "$output_dir/hist_ind_cutoffs_10.png", as(png) replace 





/*----------------------------------------------------------------------------*/

**# ******* Partial vs. NO bonus  -- 60 vs 0 *******

/*----------------------------------------------------------------------------*/

gen X=st_ind10
gen cutX=threshold_2_hg
drop if res_sned==1 // drops the 100% bonus winners



/*----------------------------------------------------------------------------*/

**# Figure 2 (a) + (b) : Histograms of the assignment variable

/*----------------------------------------------------------------------------*/



* Plots of the assignment variable
* Zeros inluded


twoway (hist X, w(0.2) fcolor(black) lcolor(black) ), ///
xtitle("Centered SNED index")   graphregion(color(white)) xlabel(-35(5)5) 
 
graph export "$output_dir/hist_ind_10_with0.png", as(png) replace 


* Plots of the assignment variable
* Zeros excluded

twoway (hist X if X!=0, w(0.2) fcolor(black) lcolor(black) ), ///
xtitle("Centered SNED index")   graphregion(color(white)) xlabel(-35(5)5)  

graph export "$output_dir/hist_ind_10_without0.png", as(png) replace 



/*----------------------------------------------------------------------------*/	

**#  Table 3(col1+2) + Fig 3 (a+b): McCrary (2008) test        

/*   https://eml.berkeley.edu/~jmccrary/DCdensity/							  */
/*----------------------------------------------------------------------------*/	
/* 

Every 3 rows for each column:
	beta= gap in the estimated distribution
	standard error of the beta
	p-value

H_0: no gap,  so looking for p-v>0.1
	
	1. 3-rows for 2006
	2. 3-rows for 2008
	3. 3-rows for 2010
	4. 3-rows for 2012	
	5. 3-rows for 2014	
	6. 3-rows for panel
*/

cd "$output_dir"

matrix define R = J(18, 2, .) // rows, cols.

matrix list R


* 2006

DCdensity X if sned_yr==2006, breakpoint(0) generate(Xj Yj r0 fhat se_fhat)   //graphname(mccrary_10_2006_with0.png)

return list

local c=1
local r=0

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))

drop Yj Xj r0 fhat se_fhat
matrix list R

 
DCdensity X if sned_yr==2006 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2006_without0.png)  

local c=2
local r=0

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* 2008

DCdensity X if sned_yr==2008, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2008_with0.png)  

return list

local c=1
local r=3

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))

drop Yj Xj r0 fhat se_fhat
matrix list R

 
DCdensity X if sned_yr==2008 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2008_without0.png)  

local c=2
local r=3

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R

 
* 2010

DCdensity X if sned_yr==2010, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2010_with0.png)  

return list

local c=1
local r=6

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))

drop Yj Xj r0 fhat se_fhat
matrix list R

 
DCdensity X if sned_yr==2010 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2010_without0.png)  

local c=2
local r=6

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* 2012

DCdensity X if  sned_yr==2012, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2012_with0.png)  
 
local c=1
local r=9

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


DCdensity X if  sned_yr==2012 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2012_without0.png)  
 
local c=2
local r=9

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* 2014

DCdensity X if  sned_yr==2014, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2014_with0.png)  
 
local c=1
local r=12

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


DCdensity X if  sned_yr==2014 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_10_2014_without0.png)  
 
local c=2
local r=12

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* panel 2006-2014


DCdensity X , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) graphname(mccrary_10_panel_with0.png)  
 
local c=1
local r=15

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


DCdensity X if X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) graphname(mccrary_10_panel_without0.png)  
 
local c=2
local r=15

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat

matrix list R

matrix colnames R = With_zeros Without_zeros

matrix rownames R = 2006-theta 2006-se 2006-pv /// 
					2008-theta 2008-se 2008-pv /// 
					2010-theta 2010-se 2010-pv /// 
					2012-theta 2012-se 2012-pv /// 
					2014-theta 2014-se 2014-pv /// 
					panel-theta panel-se panel-pv  

preserve
	clear
	svmat R
	save "$output_dir/Table3-McCrary-10.dta", replace
 
restore


/*---------------------------------------------------------------------------*/

**# Figure 1 (a) : Plots of the probability of winning against the centered assignment variable

/*---------------------------------------------------------------------------*/


** organize into bins
set more off
gen bins_X=.
sum X
forvalues i=0(0.1)25{
qui replace bins_X=`i'+0.1 if X>=`i' & X<0.1+`i' & X!=.
}

forvalues i=-0.1(-0.1)-35{
qui replace bins_X=(`i'-0.1) if X<`i' & X>=(-0.1+`i') & X!=.
}

gen 	sel35=1 if res_sned==2
replace sel35=0 if res_sned==3

collapse (mean) sel35 , by(bins_X)

twoway(scatter sel35 bins_X , msymbol(circle_hollow) mc(black) xline(0, lcolor(black))), graphregion(color(white)) /// 
xtitle("SNED index relative to homogenous group cutoffs")    graphregion(color(white)) ytitle("Pr(win SNED award)") xlabel(-35(5)5)

graph export  "$output_dir/prob_win_10.png", as(png) replace



/*----------------------------------------------------------------------------*/

**# ******* Full vs. Partial bonus -- 100 vs 60 *******

/*----------------------------------------------------------------------------*/
 
global replica_dir="<add path to your working directory>"

global data_dir=	"$replica_dir/data"
global output_dir=	"$replica_dir/output"

use "$data_dir/slvl_fin.dta", clear  



gen X=st_ind25
gen cutX=threshold_1_hg
drop if res_sned==3 // drops no winners


/*----------------------------------------------------------------------------*/

**# Figure 2 (c) + (d) : Histograms of the assignment variable

/*----------------------------------------------------------------------------*/


* Plots of the assignment variable
* Zeros included

twoway (hist X, w(0.2) fcolor(black) lcolor(black) ), ///
xtitle("Centered SNED index")   graphregion(color(white))  xlabel(-10(5)25) 

graph export "$output_dir/hist_ind_25_with0.png", as(png) replace 


* Plots of the assignment variable
* Zeros excluded

twoway (hist X if X!=0, w(0.2) fcolor(black) lcolor(black) ), ///
xtitle("Centered SNED index")   graphregion(color(white)) xlabel(-10(5)25)
  
graph export "$output_dir/hist_ind_25_without0.png", as(png) replace 



/*----------------------------------------------------------------------------*/	

**#  Table 3 (col3+4) + Fig 3 (c+d): McCrary (2008) test        

/*   https://eml.berkeley.edu/~jmccrary/DCdensity/							  */
/*----------------------------------------------------------------------------*/	
/* 

Every 3 rows for each column:
	beta= gap in the estimated distribution
	standard error of the beta
	p-value

H_0: no gap,  so looking for p-v>0.1
	
	1. 3-rows for 2006
	2. 3-rows for 2008
	3. 3-rows for 2010
	4. 3-rows for 2012	
	5. 3-rows for 2014	
	6. 3-rows for panel
*/

cd "$output_dir"

matrix define R = J(18, 2, .) // rows, cols.

matrix list R


* 2006

DCdensity X if sned_yr==2006, breakpoint(0) generate(Xj Yj r0 fhat se_fhat)  // graphname(mccrary_25_2006_with0.png)   

return list

local c=1
local r=0

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))

drop Yj Xj r0 fhat se_fhat
matrix list R

 
DCdensity X if sned_yr==2006 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2006_without0.png)  

local c=2
local r=0

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R



* 2008

DCdensity X if sned_yr==2008, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2008_with0.png)  

return list

local c=1
local r=3

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))

drop Yj Xj r0 fhat se_fhat
matrix list R

 
DCdensity X if sned_yr==2008 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2008_without0.png)  

local c=2
local r=3

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R

 
* 2010

DCdensity X if sned_yr==2010, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2010_with0.png)  

return list

local c=1
local r=6

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))

drop Yj Xj r0 fhat se_fhat
matrix list R

 
DCdensity X if sned_yr==2010 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2010_without0.png)  

local c=2
local r=6

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* 2012

DCdensity X if  sned_yr==2012, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2012_with0.png)  
 
local c=1
local r=9

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


DCdensity X if  sned_yr==2012 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2012_without0.png)  
 
local c=2
local r=9

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* 2014

DCdensity X if  sned_yr==2014, breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2014_with0.png)  
 
local c=1
local r=12

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


DCdensity X if  sned_yr==2014 & X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) //graphname(mccrary_25_2014_without0.png)  
 
local c=2
local r=12

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


* panel 2006-2014


DCdensity X , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) graphname(mccrary_25_panel_with0.png)  
 
local c=1
local r=15

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat
matrix list R


DCdensity X if X!=0 , breakpoint(0) generate(Xj Yj r0 fhat se_fhat) graphname(mccrary_25_panel_without0.png)  
 
local c=2
local r=15

matrix R[`r'+1, `c']= r(theta)
matrix R[`r'+2, `c']= r(se)
matrix R[`r'+3, `c'] = 2 * normal(-abs(R[`r'+1, `c'] / R[`r'+2, `c']))
drop Yj Xj r0 fhat se_fhat

matrix list R

matrix colnames R = With_zeros Without_zeros

matrix rownames R = 2006-theta 2006-se 2006-pv /// 
					2008-theta 2008-se 2008-pv /// 
					2010-theta 2010-se 2010-pv /// 
					2012-theta 2012-se 2012-pv /// 
					2014-theta 2014-se 2014-pv /// 
					panel-theta panel-se panel-pv  

preserve
	clear
	svmat R
	save "$output_dir/Table3-McCrary-25.dta", replace
 
restore


/*---------------------------------------------------------------------------*/

**# Figure 1 (b) : Plots of the probability of winning against the centered assignment variable

/*---------------------------------------------------------------------------*/


** organize into bins
set more off
gen bins_X=.
sum X
forvalues i=0(0.1)25{
qui replace bins_X=`i'+0.1 if X>=`i' & X<0.1+`i' & X!=.
}

forvalues i=-0.1(-0.1)-35{
qui replace bins_X=(`i'-0.1) if X<`i' & X>=(-0.1+`i') & X!=.
}

gen 	sel35=1 if res_sned==1
replace sel35=0 if res_sned==2

collapse (mean) sel35 , by(bins_X)

twoway(scatter sel35 bins_X , msymbol(circle_hollow) mc(black) xline(0, lcolor(black))), graphregion(color(white)) /// 
xtitle("SNED index relative to homogenous group cutoffs")    graphregion(color(white)) ytitle("Pr(win SNED award)") xlabel(-10(5)20) 

graph export "$output_dir/prob_win_25.png", as(png) replace


*----------------------------------------------------------------------------END



